#define _CRT_SECURE_NO_DEPRECATE
#include <cstdio>
#include <fcntl.h> // for _O_BINARY
#include <vector>
#include <string>
#include <cstring>
#include <algorithm>
#include <cmath>
#include "grate_constants.h"
#include "ext_rational.h"
#include "words_queue.h"
#include "dfa.h"
#include "dfa_construct.h"
#include "ad_construct.h"
#include "input_data.h"
#include "grate_config.h"

using namespace std;

// Contains all the parameters
grate_config CONFIG;

// Main containers - queue and DFA table
words_queue * WORDS_QUEUE;
dfa * DFA;

//------------------------------------------- INITIALIZATION -----------------------------------

// Determines if queue size is more important than dfa size
int queue_more_important(input_data & INPUT){
  int alph_size = INPUT.ALPH_SIZE;
  if (alph_size <= 2) return 0;
  return (INPUT.beta < ext_rational(alph_size-1, alph_size-2, 1));
}

void create_queue(input_data & INPUT, int size){
  int qsize;
  if (size != -1){
    qsize = size;
  }
  else{ // define size automatically
    double max_size = INPUT.ROOT_LENGTH * pow((double)INPUT.ALPH_SIZE, (double)INPUT.ROOT_LENGTH);
    if (sizeof(char) * max_size < LITTLE_MEMORY){
      qsize = (int)max_size + 1;
    }
    else{
      qsize = QUEUE_MEMORY_STEP;
      while (1){
        try{
          char * buff = new char [qsize + QUEUE_MEMORY_STEP];
          delete[] buff;
          qsize += QUEUE_MEMORY_STEP;
        }
        catch(...){
          break;
        }
      }
      qsize >>= 1;
    }
  }
  bool use_long_mask = (INPUT.ALPH_SIZE > 4);
  WORDS_QUEUE = new words_queue(qsize, use_long_mask);
}

void create_dfa(input_data & INPUT, int size){
  int dsize;
  if (size != -1){
    dsize = size;
  }
  else{ // define size automatically
    double max_size = INPUT.ROOT_LENGTH * pow((double)INPUT.ALPH_SIZE, (double)INPUT.ROOT_LENGTH);
    if (sizeof(int) * max_size * INPUT.ALPH_SIZE < LITTLE_MEMORY){
      dsize = (int)max_size + 1;
    }
    else{
      dsize = DFA_MEMORY_STEP;
      while (1){
        try{
          int * buff = new int [(dsize + DFA_MEMORY_STEP) * INPUT.ALPH_SIZE];
          delete[] buff;
          dsize += DFA_MEMORY_STEP;
        }
        catch(...){
          break;
        }
      }
      dsize >>= 1;
    }
  }
  DFA = new dfa(dsize, INPUT.ALPH_SIZE);
}

void init(input_data & INPUT){
  if (CONFIG.auto_queue_size == false){
    create_queue(INPUT, CONFIG.QUEUE_SIZE);
  }
  if (CONFIG.auto_dfa_size == false){
    create_dfa(INPUT, CONFIG.DFA_SIZE);
  }
  
  if (queue_more_important(INPUT)){
    if (CONFIG.auto_queue_size) create_queue(INPUT, -1);
    if (CONFIG.auto_dfa_size) create_dfa(INPUT, -1);
  }
  else{
    if (CONFIG.auto_dfa_size) create_dfa(INPUT, -1);
    if (CONFIG.auto_queue_size) create_queue(INPUT, -1);
  }

  if (CONFIG.OUTPUT_PROGRESS){
    if (CONFIG.auto_queue_size) INPUT.dump("Amount of memory allocated for queue = %d bytes\n", sizeof(char)*WORDS_QUEUE->size);
    if (CONFIG.auto_dfa_size) INPUT.dump("Amount of memory allocated for DFA = %d bytes\n", sizeof(int)*(DFA->size)*INPUT.ALPH_SIZE);
  }
}


//--------------------------------------------------- ANSWER CALCULATION-----------------------------------------

double * cnt_new;
double * cnt_old;

double grate;

void cnt_algo2(input_data & INPUT){
  int i, j, k;
  char h;
  int total = DFA->total;
  cnt_old = new double [total+1];
  memset(cnt_old, 0, sizeof(double) * (total + 1));
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("old counters allocated\n");
  }
  cnt_new = new double [total+1];
  memset(cnt_new, 0, sizeof(double) * (total + 1));
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("new counters allocated\n");
  }
  cnt_old[1] = 1;
  double pcomp = 1, ccomp, last_loc_max = 0, last_loc_min = 0;
  double * rates = new double [CONFIG.MAXITER];
  for (i=0; i<CONFIG.MAXITER; i++){
    for (j=1; j<=total; j++){
      for (h=0; h<INPUT.ALPH_SIZE; h++){
        k = (*DFA)[j][h];
        if (!k) continue;
        cnt_new[k] += cnt_old[j];
      }
    }
    ccomp = 0;
    for (j=1; j<=total; j++){
      if (CONFIG.ADD_LOOPS) cnt_new[j] += cnt_old[j];
      ccomp += cnt_new[j];
      cnt_old[j] = cnt_new[j];
      cnt_new[j] = 0;
    }
    rates[i] = ccomp / pcomp;
    // ORACLE HERE
    if (i >= MINITER){
      if ((rates[i-1] >= rates[i-2]) && (rates[i-1] >= rates[i])) last_loc_max = rates[i-1];
      if ((rates[i-1] <= rates[i-2]) && (rates[i-1] <= rates[i])) last_loc_min = rates[i-1];
      if (fabs(rates[i] - rates[i-1]) < INPUT.delta){
        // Oscillating
        if ((last_loc_max > 0) && (last_loc_min > 0)){
          if (last_loc_max - last_loc_min < INPUT.delta){
            if (CONFIG.OUTPUT_PROGRESS){
              INPUT.dump("Oscillating sequence! The number of iterations = %d\n", i);
            }
            if (GUI_VERSION){
              printf("%d\n", i);
            }
            grate = rates[i] - (CONFIG.ADD_LOOPS == true);
            return;
          }
        }
        // Monotonic
        else{
          if (CONFIG.OUTPUT_PROGRESS){
            INPUT.dump("Monotonic sequence! The number of iterations = %d\n", i);
          }
          if (GUI_VERSION){
            printf("%d\n", i);
          }
          double q = pow(fabs(rates[i]-rates[i-1]), 1.0 / i);
          double sum = (rates[i] - rates[i-1]) / (1 - q);
          grate = rates[i-1] + sum - (CONFIG.ADD_LOOPS == true);
          return;
        }
      }
    }
    pcomp = ccomp;
  }
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("\nNot enough iterations! Current difference is %.20lf\n\n", rates[i-1] - rates[i-2]);
  }
  if (GUI_VERSION) printf("-1\n");
  grate = rates[i-1] - (CONFIG.ADD_LOOPS == true);
}

// Edge list
int * edge1;
int * edge2;
bool have_edge_list;

// Some values to output
int scc_num, nontriv_scc_num, scc_size, scc_edgnum;

int * zcc;

void cnt_algo1(input_data & INPUT, int id){
  double minrate, maxrate, crate, alpha;
  int i, iter = 0;
  while (iter < CONFIG.MAXITER){
    iter++;
    if (have_edge_list){
      for (i=0; i<scc_edgnum; i++){
        cnt_new[edge2[i]] += cnt_old[edge1[i]];
      }
    }
    else{
      /*for (i=1; i<=scc_size; i++){
        for (j=0; j<ALPH_SIZE; j++){
          if (dfa[i*ALPH_SIZE + j] > 0){
            cnt_new[dfa[i*ALPH_SIZE + j]] += cnt_old[i];
          }
        }
      }*/
    }
    minrate = 1e100;
    maxrate = 0;
    for (i=1; i<=scc_size; i++){
      if (CONFIG.ADD_LOOPS){
        cnt_new[i] += cnt_old[i];
      }
      crate = cnt_new[i] / cnt_old[i];
      if (crate < minrate) minrate = crate;
      if (crate > maxrate) maxrate = crate;
      cnt_old[i] = cnt_new[i];
      cnt_new[i] = 0;
    }
    if (maxrate - minrate < 2 * INPUT.delta){
      alpha = (minrate + maxrate) / 2 - (CONFIG.ADD_LOOPS == true);
      if (minrate > MAX_GRATE) alpha = 0;
      if (alpha > grate) grate = alpha;
      if (CONFIG.OUTPUT_PROGRESS) INPUT.dump("\n%d iterations used\n\n", iter);
      INPUT.dump("Growth rate at scc %d: %.20lf\n", id, alpha);
      if (GUI_VERSION) printf("%d\n", iter);
      return;
    }
  }
  if (CONFIG.OUTPUT_PROGRESS) INPUT.dump("\nNot enough iterations or precision error! %.20lf precision only!\n\n", (maxrate-minrate)/2);
  if (GUI_VERSION) printf("-1\n");
  alpha = (maxrate + minrate) / 2 - (CONFIG.ADD_LOOPS == true);
  INPUT.dump("Growth rate at scc %d: %.20lf\n", id, alpha);
  if (alpha > grate) grate = alpha;
  
}

void solve_scc(input_data & INPUT, int id){
  int i, j;
  int total = DFA->total;

  scc_edgnum = 0;
  scc_size = 0;
  
  for (i=1; i<=total; i++){
    if (zcc[i] == id){
      scc_size++;
      zcc[i] = -scc_size;
    }
  }

  for (i=1; i<=total; i++){
    if (zcc[i] < 0){
      for (j=0; j<INPUT.ALPH_SIZE; j++){
        int v = (*DFA)[i][j];
        if (v > 0 && zcc[v] < 0){
          scc_edgnum++;
        }
      }
    }
  }

  have_edge_list = true;
  
  edge1 = new int [scc_edgnum];
  edge2 = new int [scc_edgnum];
  int cnum = 0, znum = 0;
  for (i=1; i<=total; i++){
    if (zcc[i] < 0){
      for (j=0; j<INPUT.ALPH_SIZE; j++){
        int v = (*DFA)[i][j];
        if (v > 0 && zcc[v] < 0){
          edge1[cnum] = -zcc[i];
          edge2[cnum++] = -zcc[v];
        }
      }
    }
  }

  for (i = 1; i <= total; i++){
    if (zcc[i] < 0) zcc[i] = 0;
  }

  cnt_old = new double [scc_size + 1];
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("old counters allocated\n");
  }
  cnt_new = new double [scc_size + 1];
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("new counters allocated\n");
  }
  
  for (i=1; i<=scc_size; i++){
    cnt_new[i] = 0;
    cnt_old[i] = 1;
  }
  
  if (scc_edgnum == scc_size){
    if (grate == 0) grate = 1;
    if (have_edge_list){
      delete[] edge1;
      delete[] edge2;
    }
    return;
  }
  
  nontriv_scc_num++;
  cnt_algo1(INPUT, id);
  if (have_edge_list){
    delete[] edge1;
    delete[] edge2;
  }
}

// TARJAN ALGORITHM

int * mtime;
int * father;
int * stack;
int ssize;
char * nxt_used; // binary encodeing here: 15-mask of number of edges done [+16] = used
char * scc;

void dfs(input_data & INPUT){
  int total = DFA->total;
  memset(nxt_used, 0, sizeof(char) * (total + 1));
  memset(mtime, 0, sizeof(int) * (total + 1));
  ssize = 0;
  father[1] = 0;
  mtime[1] = 1;
  int ctime = 1, cver = 1, cedge, nver, pver;
  while (cver > 0){
    cedge = nxt_used[cver] & LONG_MASK;
    if (cedge == 0){ 
      nxt_used[cver] ^= 16; // is in stack
      stack[++ssize] = cver;
    }
    else{ // check the son
      pver = (*DFA)[cver][cedge-1];
      if ((pver > 0) && (father[pver] == cver) && (mtime[pver] < mtime[cver])){
        mtime[cver] = mtime[pver];
        nxt_used[cver] |= 32; // improvement
      }
    }
    if (cedge >= INPUT.ALPH_SIZE){
      // handle vertex away

      if ((nxt_used[cver] & 32) == 0){ // no improvement
        if (stack[ssize] == cver){
          nxt_used[cver] ^= 16; // ouf of stack
          ssize--;
        }
        else{
          scc_num++;
          do{
            pver = stack[ssize--];
            nxt_used[pver] ^= 64; // in scc
            scc_size++;
            nxt_used[pver] ^= 16; // ouf of stack
            scc[pver] = scc_num;
          }while (pver != cver);
        }
      }

      cver = father[cver];
      continue;
    }
    nxt_used[cver]++;
    nver = (*DFA)[cver][cedge];
    if (!nver) continue;
    if (!mtime[nver]){
      ctime++;
      mtime[nver] = ctime;
      father[nver] = cver;
      cver = nver;
      continue;
    }
    else if ((nxt_used[nver] & 16) && (mtime[nver] < mtime[cver])){
      mtime[cver] = mtime[nver];
      nxt_used[cver] |= 32; // improvement
    }
  }
}

void divide_scc_solve(input_data & INPUT){
  int total = DFA->total;
  mtime = new int [total + 1];
  father = new int [total + 1];
  stack = new int [total + 1];
  nxt_used = new char [total + 1];
  
  scc = new char [total + 1];
  memset(scc, 0, sizeof(char) * (total + 1));
  
  dfs(INPUT);
  
  delete[] mtime;
  delete[] stack;
  delete[] father;
  delete[] nxt_used;
  
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("The number of non-singleton strongly connected components: %d\n", scc_num);
  }

  zcc = new int[total + 1];
  for (int i = 1; i <= total; i++){
    zcc[i] = scc[i];
  }
  delete[] scc;

  for (int cur_scc=1; cur_scc<=scc_num; cur_scc++){
    solve_scc(INPUT, cur_scc);
    if (CONFIG.OUTPUT_PROGRESS){
      INPUT.dump("The number of vertices in scc %d: %d\n", cur_scc, scc_size);
      INPUT.dump("The number of edges in scc %d: %d\n", cur_scc, scc_edgnum);
      INPUT.dump("\n");
    }
  }

  delete[] zcc;

  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("The number of non-trivial scc: %d\n", nontriv_scc_num);
  }
}

// ------------------------------------ CLEAR DATA -------------------------

void clear_data(){
  grate = 0;
  scc_num = 0;
  nontriv_scc_num = 0;
  scc_size = 0;
  scc_edgnum = 0;
  
  /*if (USE_ORACLE){
    delete[] dfa;
  }
  if (!USE_ORACLE && !have_edge_list){
    delete[] dfa;
  }*/
  delete DFA;

  delete[] cnt_new;
  delete[] cnt_old;
}

void solve(input_data & INPUT){

  INPUT.create_logfile();
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump_input();
  }

  init(INPUT);

  //build_antidictionary();
  i64 max_words_queue_size = construct_antidictionary(INPUT, WORDS_QUEUE, DFA);

  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("Maximal queue size = %I64d\n", max_words_queue_size);
    INPUT.dump("The number of words = %d\n", DFA->words_num);
    INPUT.dump("The number of vertices = %d\n", DFA->total);
  }

  /*if (GUI_VERSION){
    if (use_long_mask) max_words_queue_size >>= 1;
    else max_words_queue_size >>= 2;
    printf("%I64d\n", max_words_queue_size);
    printf("%d\n", words_num);
    printf("%d\n", total);
  }*/

  construct_dfa(DFA);

  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("\nDFA built!\n\n");
  }

  if (CONFIG.USE_ORACLE){
    cnt_algo2(INPUT);
  }
  else{
    divide_scc_solve(INPUT);
  }

  INPUT.dump("Upper bound: %.13lf\n", grate);
  printf("%.13lf\n", grate);

  fflush(stdout);
  fclose(INPUT.logfile);
  clear_data();
}

int main(int argc, char * argv[]){
  CONFIG.parse_config_file();
  input_data INPUT;
  // File IO
  if (argc >= 2){
    freopen(argv[1], "r", stdin);
    if (argc >= 3){
      freopen(argv[2], "w", stdout);
    }
    while (scanf("%d%d", &INPUT.ALPH_SIZE, &INPUT.ROOT_LENGTH) > 0){
      if (INPUT.ALPH_SIZE < 0) break;
      scanf("%d%d%d", &INPUT.beta.num, &INPUT.beta.denum, &INPUT.beta.p);
      scanf("%d", &INPUT.power_type);
      scanf("%lf", &INPUT.delta);
      solve(INPUT);
    }
  }
  // Or console mode
  else{
    printf("Hello! Working in console mode now...\n");
    while (1){
      printf("Alphabet size: ");
      scanf("%d", &INPUT.ALPH_SIZE);
      printf("Max length of the root: ");
      scanf("%d", &INPUT.ROOT_LENGTH);
      printf("Beta numerator: ");
      scanf("%d", &INPUT.beta.num);
      printf("Beta denominator: ");
      scanf("%d", &INPUT.beta.denum);
      printf("Is there a + ?: ");
      scanf("%d", &INPUT.beta.p);
      printf("Type of abelian power (0 - weak, 1 - moderate, 2 - strong: ");
      scanf("%d", &INPUT.power_type);
      printf("Delta (exponential form is acceptable): ");
      scanf("%lf", &INPUT.delta);
      solve(INPUT);
      printf("\nContinue? (0/1): ");
      int ch;
      scanf("%d", &ch);
      if (!ch) break;
    }
  }
  return 0;
}
